Quantitative modeling identifies critical cell mechanics driving bile duct lumen formation

Biliary ducts collect bile from liver lobules, the smallest functional and anatomical units of liver, and carry it to the gallbladder. Disruptions in this process caused by defective embryonic development, or through ductal reaction in liver disease have a major impact on life quality and survival of patients. A deep understanding of the processes underlying bile duct lumen formation is crucial to identify intervention points to avoid or treat the appearance of defective bile ducts. Several hypotheses have been proposed to characterize the biophysical mechanisms driving initial bile duct lumen formation during embryogenesis. Here, guided by the quantification of morphological features and expression of genes in bile ducts from embryonic mouse liver, we sharpened these hypotheses and collected data to develop a high resolution individual cell-based computational model that enables to test alternative hypotheses in silico. This model permits realistic simulations of tissue and cell mechanics at sub-cellular scale. Our simulations suggest that successful bile duct lumen formation requires a simultaneous contribution of directed cell division of cholangiocytes, local osmotic effects generated by salt excretion in the lumen, and temporally-controlled differentiation of hepatoblasts to cholangiocytes, with apical constriction of cholangiocytes only moderately affecting luminal size.

where γ denotes the friction coefficient, l 0 ij = ||r 0 ij || = ||r 0 j − r 0 i || and l ij = ||r ij || are the initial (cell at rest) and actual lengths between the nodes, and v ij = v j − v i is the relative velocity of nodes i, j. The force balance equation with external forces F ext demands that F ext + F int = 0, hence: The linear spring constant for a sixfold symmetric triangulated lattice can be related approximately to the cortex Young modulus E cor with thickness h cor by [2] Besides tension, the cortex resists to bending. The bending resistance in the cortex is incorporated by the angular resistance of the hinges determined by two adjacent triangles T 1 = {ijk} and T 2 = {ijl}. This can be accounted for by the bending moment M : where k b is the bending constant, and θ is is the angle between the normal vectors to the triangles n α , n β by their scalar product (n α n β ) = cos(θ). θ 0 is the angle of spontaneous curvature. A spontaneous curvature denotes the curvature for which the bending energy is zero. The moment M can be transformed to an equivalent force system F m,z (z ∈ {ijl}) for the triangles T 1 and T 2 where here for T 1 we can compute F m,i = M/l 1 n 1 using l 1 as the distance between the hinge of the triangle pair and the point i, and similar expression for F m,l . The forces working on nodes j, k must at least fulfill F m,i + F m,j + F m,k + F m,l = 0 to conserve total momentum, see e.g. [3,4]. The bending stiffness of the cortex is approximated by where ν ≈ 0.5 is the Poisson's ratio of the cortex. Note that the elastic modulus of the cortex, E cor enters both, the bending force and the tension force.

Volume forces
The volume change of the cell depends on the applied pressure and the cell bulk modulus K V . The total compressibility of the cell depends on volume fraction of water in the cytosol, the cytoskeleton (CSK) volume fraction and structure, and the compressibility of the individual organelles. In addition, it may be influenced by the permeability of the plasma membrane for water, the presence of caveolae, and active responses inside and of the cell. We calculate the internal pressure in a cell by the logarithmic strain for volume change: whereby V is the actual volume and V 0 is the reference volume i.e., the volume of the cell not subject to compression forces. For small deviations of V from V 0 , Within our model the volume V of the cell is computed summing up the volumes of the individual tetrahedra that build up the cell. The nodal force is obtained by multiplying the pressure with the nodal Voronoi area S i (see [5]), i.e. F vol,i = pS i R i where the local curvature vector R i is computed for that node using a discrete Laplace-Beltrami operator [7].

Contact model and adhesive forces
Whereas in a CBM, cells interact through central forces described by (modified) Hertz or JKR theory for adhesive spheres, in DCM the interaction forces need to be defined for each node individually, thereby requiring a representation of local surface heterogeneities. The approach followed in this work adopts a discretised Maugis-Dugdale theory. The Maugis-Dugdale theory for adhering bodies is a generalization of the JKR theory for spheres [6]. This theory captures the full range between the Derjaguin-Muller-Toporov (DMT) zone of long reaching adhesive forces of a hard homogeneous isotropic elastic sphere and small adhesive deformations in the Johnson-Kendall-Roberts (JKR)-limit of a soft homogeneous isotropic elastic sphere of short interaction ranges. Here, we associate to each triangle of the cell surface a circumscribing sphere reflecting the local curvature. Two triangles belonging to different cells can locally interact by collision of their associated spheres. To compute the magnitude of these interactions with respect to the nodes of the triangles, we use the general Maugis-Dugdale stress, which is integrated over the common contact area between the triangles, resulting in the net forces on the nodes of the two interacting triangles (see [3,7] for more details). The adhesion force is fully determined by the specific adhesion energy W and the typical effective adhesive range h 0 that reflects the attractive cutoff distance between the bodies. We set h 0 = 2 · 10 −8 m in all the simulations [8].
In our model, Maugis-Dugdale theory is applied to every set of triangles which constitute the cell surface. A varying number of cadherin bonds is mimicked by varying the adhesion energy along the cell surface triangles.

Migration forces
In models in which cell movement and deformation is mimicked by force-balance equations following Newton's law of motion as this is the case for the DCM and CBM, migration is usually modelled by an active migration force F mig representing the random micro-motility of a cell. For the sake of simplicity the cell shape, filopodia etc. movements during migration are not detailed in the model but instead the different mechanisms in cell migration is lumped into one net force which is uniformly distributed to the nodes the cell if not mentioned otherwise. For specific applications, the forces might be non-uniformly distributed. In absence of influences that impose a certain direction or persistence, it is commonly assumed that the migration force is stochastic, formally resulting in F mig = F ran , with F ran = 0, and F ran (t) ⊗ F ran (t ) = Mδ(t − t ), where M is an amplitude 3 × 3 matrix and relates to the diffusion tensor D of the cell. As cell migration is active, depending on the local matrix density and orientation of matrix fibers, the autocorrelation amplitude matrix M cannot a priori be assumed to follow a fluctuation-dissipation (FD) theorem. However, "measuring" the position of a cell in the simulations the position autocorrelation function might be experimentally used to determine the diffusion tensor using ((r(t + τ ) − r(t)) ⊗ (r(t + τ ) − r(t)) = 6Dτ , and M be calibrated such that the numerical solution of the equation of motion for the cell position reproduces the experimental result for the position-autocorrelation function. For example, in a homogeneous environment M can be casted into a form formally equivalent to the FD-theorem, leading to an k B T -equivalent for cellular systems, that is controlled by the cell itself [9,10].
Note, we assume here a momentum transfer to the ECM by applying the ECM friction and active micro-motility term, but we do not model the ECM explicitly.

Cell cycle and cell division
During cytokinesis, the continuous shrinking of the contractile ring, together with the separation of the mitotic spindle, gradually divides the mother cell into two daughter cells. After mitosis the cell has split up in two adhering cells. Such a process of cell division in 2D deformable cells has been previously described (e.g. [11]) but it can be costly and tedious to perform in 3D.
As we are merely interested in long term effects (i.e. hundreds of cell divisions), and as the cytokinesis is a short process compared to the duration of the entire cell cycle, we avoid these particular tedious intermediate steps in our model, and directly create two new adhering cells that are within the shape of the mother cells just before its division using the basic algorithm established in ref. [12].
In our algorithm, the two created daughter cells are initially round and separated from each other by the division plane. They grow artificially fast 1 until they reach the boundary of the mother envelope. The division plane is removed from the system after the daughter cells have found a mechanical equilibrium in this envelope. We have here slightly extended our algorithm to mimic the very short rounding-up period just before division (see Fig 2G in main text). This is achieved by setting the equilibrium edge lengths (see below) of the triangulation temporarily to a smaller value. This can cause a slight increased pressure in the cell during division.
In the simulations, cells grow by increasing their volume and surface and they can divide when their actual volume reaches a target volume. This target volume is predefined and is for a newborn cell by default set to twice the volume of that newborn cell. The orientation of cell division can be chosen randomly (as would be the case for non polarized hepatocytes) or with a preferred direction. The latter is the case when a polarized cell divides. In this case the division direction is the polarization vector. As observed experimentally, we also implement a temporary and short rounding up period just before cells division. This has a consequence that the cell volume slightly drops and hence the pressure inside the cell increases, causing a rounding of the surface. The shortening of the vertex edges is chosen such that this pressure increase is not more than 100 to 200 P a, achieved by modifying the lengths from l 0 to 0.95 * l 0 . The rounding up period in our model does not take more than 5% of the total cell cycle time.
The cells in our system can either be proliferating or quiescent. Immediately after cell division, all cells are assumed to proliferate. A change to quiescence is here assumed to be induced by chemical signals of other cells. To ensure in the simulations that only the fraction of the cells is growing conform with the experimental observations, an algorithm is invoked that sweeps every time step over a certain cell type and ensures that the prescribed fraction of proliferating cells is maintained. The picking of cells that need to go to quiescence is done randomly. This algorithm was chosen as the mechanistic control of how many cells divide was not subject of the present simulation study.

DCM extension: Polarity vector
A normalized polarity vector P (PCP) is assigned to each of the cells. This vector defines for every cell two opposite conical polar regions (see Fig 2A in main text ). The triangular regions on the cell surface that are marked as polar are defined by the scalar product condition: where N i is a normalized vector with origin the center of mass of the cell, and pointing to the middle of triangle i. The area of the polar region a can be chosen according to 0 < a < 1. We chose in the simulations of this work a = 0.7. This corresponds to a bi-conical area with angle α of approximately 60 degrees. The polar regions (triangles) may have different physical properties than the rest of the cell surface, such as a different specific adhesion energy. The PCP vector is assumed to be perpendicular to the apico-basal vector, which is enforced in the model with every timestep. (We note here that in our planar simulations neither the PCP or ABP have a Z component.)

DCM extensions : apical vector, apical constriction and tight junctions
Every cell also has an apical vector A (ABP), which defines the apical and basal side of the cell. Similar to the polar regions, the apical and basal triangular regions on the cell surface are defined by the conditions: and Like for the polar vector, this corresponds to a conical area with angle β and the scalar value of b ∈ [0, 1] determines the area of the regions, however the difference is that the apical side of the cell may have different physical properties than the basal side. In this work the cortical cytoskeleton on the apical side can contract, causing apical constriction.
In the model, the apical vector is assumed every time step to point towards the principal region where the cell surface has the most free area, in line with the assumption made that the apical side forms the cell's interface with the lumen. Because for every triangle of a cell the free area and contact area is known after applying the contact model, the principal vector A can obtained by Here the summation only counts triangles that have a very small contact area with a triangle from another cell. These conditions can for example be fulfilled when gaps between cells are created due to osmotic effects (e.g. onset of a small cavity). A is updated every timestep, except during cell division and a short relaxation phase right after cell division.
Adherent cells can also develop tight junctions (TJ) between them. In the model, we define for every cell that has an apical vector, a region with TJ triangles i that fulfill the condition : Here c 1 and c 2 are the limit values that define a conical ribbon along the cell where TJ can be present (see Fig 2D in main text). The parameter values for c 1 and c 2 are chosen such that the ribbon is about 1 or maximum at 2 triangles thick. As TJs denote areas on the cell surface with a reinforced connection, a larger specific adhesion energy is assigned to these areas than to the other surface areas. As such, when two triangles of different cells with a TJ come into contact, the force necessary to separate them will be much higher than for normal cell-cell contacts.
The apical constriction on the apical side is modeled as a shortening of the equilibrium element lengths between connected nodes that are located within the apical region (see below). This results in a movement of the two nodes towards each other. The shortening operation is executed only once. After a certain relaxation time a new force equilibrium is reached between the nodes. The modification of the equilibrium element lengths of the nodal springs during apical constriction are computed assuming two different zones: which is the zone of the tight junctions, and in the apical domain.
To simulate the presence of an apical circumferential ring, the strongest constriction is applied in the zone of the TJ ribbon, whereas a lower constriction is applied in the apical domain (0 < d 2 < d 1 < 1). We have adopted the values d 2 = 0.5 and d 1 = 0.7. This warrants a strong contraction 2 . Moreover, the adhesive energy on the apical side is assumed to be lower than elsewhere on the cell because generally one observes that cadherin staining is less strong there (see Fig A1B).

DCM extensions : Osmotic effects and signalling
We assume that osmotic forces in the system result from differences of solutes e.g by molecules such as salts that are released from the cell into the extracellular space. These molecules further diffuse locally and cause a concentration gradient with more Figure A1. Cartoon indicating the parameters needed to determine the apical domain and apical circumferential ring. distant locations. For example, if one cell starts to excrete ions, a salt concentration gradient is generated attracting water molecules towards the cell potentially causing a hydrostatic pressure increase strong enough to push away cells in the direct neighborhood of the central cell. As a consequence an extracellular space can arise, depending in shape and size on the dynamics and mechanics of the surrounding cells. One way to simulate the diffusion process of salts and flow of water could be to solve a system of partial differential equations (PDEs) for diffusive transport of salt and the advection of water. As the surrounding cells reorganise during the lumen formation process and hence the shape of the extracellular space changes, the boundary conditions for the PDEs change constantly. The constant change requires high resolution meshing and constant re-meshing of the complex cavity, which is numerically complicated and computationally time consuming. To alleviate this, inspired by smoothed particle hydrodynamics methods [13], we introduce the concept of "tracer particles (TP)", which are small "inert" particles. They can diffuse freely inside the cell or in the extracellular space without mutual interaction (see Fig 2F in main text). All the tracer particles (we used about 500 in numbers per cell) are initially inside the cell and start diffusing from there. The advantage of using particles to mimic the fluid is that their motion can be simulated by a Newton's equation of motion as for the cells, hence the coupling of fluid movement and cell deformation is straightforward. The force on each particle originates principally from Brownian motion as and its magnitude is controlled by the TP diffusion coefficient. Additionally, a small force component is exerted to them in the apical direction to facilitate their movement towards the apical side which can be justified by directional vesicle transport towards lumen in polarized cells [14].
To ensure that the TP can leave the cell only through a certain region, we define a "transparent" region for each cell surface through which these particles can move into the extracellular space. The transparent region is equal to the apical side of the cells. In other parts of the cell, the particles cannot move across the cell surface. The particles can thus locally cross the cell boundary when moving from inside to outside, but are always repelled when trying to move back from the outside to the inside of the cell. These particles do not represent individual ions but rather each population of ions approximately obtained by integrating the spatial density of ions over local space regions. Nevertheless, to ensure diffusion constant is similar to that of the ions, the particles are assigned with the diffusion coefficient of the ion species.
In the extracellular space a tracer particles marks a triangle if it is in close proximity of that triangle representing a free cell surface piece, not in full contact with another cell. As there are many TP present, all free surface areas of a cell will progressively become marked. To simulate the resulting hydrostatic pressure on that cell, each marked triangle receives an external force F osm = P A i n i , where P is the assumed hydrostatic pressure that corresponds to the osmotic pressure generated by the difference in ion concentration (assumed constant here), and A i is the surface of the marked triangle with local normal n i . Triangles belonging to two different cells that are in contact with each other have a negligible probability to be marked as the particles cannot access the space between them.

CBM forces
The CBM does not resolve cell shape explicitly. Forces are in the CBM are assumed to be exerted on the center of mass of the cell.

Adhesive and repulsive forces
In the CBM, cells are approximated by homogeneous, isotropic, elastic and adhesive spheres which split into two adherent cells during mitosis. Under conditions met in this work [16,17], the total cell to cell interaction force can be approximated by the JKR force. The interaction force is computed by The contact radius a in Eq. 14 if a function of the overlap δ ij = || r j − r i || − R i − R j between the cell, allows to compute the cell-cell contact area, and can be obtained from: In the latter equations,Ê andR are defined aŝ with E i and E j being the cell Young's moduli, ν i and ν j the Poisson numbers and R i and R j the radii of the cells i and j, respectively. To enforce consistency with the by construction more accurate DCM that explicitly accounts for multi-body interactions, the Young modulus for every cell E i in the JKR model was replaced by an "apparent" Young's modulusẼ i that increases as function of the local cell density [12]: The parameters a i can be determined in a calibration experiment, where one compresses a spheroid and measures the force between the cells [1,12].
Because the CBM and the DCM are used in hybrid mode, they can make contact with each other. Using the Maugis-Dugdale formulation this can be solved in a natural way. The contact between a CBM and a triangle of the DCM is resolved using the triangle-triangle contact formulation for two DCMs, yet one of the triangles circumscribing sphere is here replaced by the sphere of the CBM cell itself. When the contact force is calculated, it is for the CBM directly passed to its center of mass.

CBM: Migration forces
The same principles for modeling the migration force in the DCM are applied to the CBM, whereby the migration force in the CBM is directly applied to the cell center.

CBM: Cell division
As in the DCM, if the cell passed a critical volume V crit = 2V 0,i , the cell undergoes mitosis and two new cells with volume V 0,i /2 are created. A simplified version of the division algorithm consists of placing directly two smaller daughter cells in the space originally filled by the mother cell at the end of the interphase [15,16]. When the two daughter cells are created, their reorganize in space driven by their cell-cell interaction force as well as the interaction forces of the two daughter cells with their other surrounding cells until mechanical equilibrium is reached. If the space filled by the mother cell is small, which is often the case for cells in the interior of a cell population, the local interaction forces occurring after replacing the mother cell by two spherical daughters, can adopt large (un-physiological) values leading to unrealistic large cell displacements. This can be circumvented by truncating the contact force to a maximal value, or temporarily reducing the contact stiffness between the daughter cells during division.

Cell-to-cell signalling
In our model, we assume that Notch-jagged signaling between cholangiocytes and hepatocytes is governed by several of conditions. First, we assume that their common contact area must be sufficiently high. Thus, a minimal number of triangles (or surface area) of both cells must be in mutual contact: A T > A min , where A T is the total contact area between two cells, and the parameter A min is set arbitrarily to 20% of the cell surface. We found that this parameter does not influence the results significantly provided that A min is larger than zero. However, this condition alone is not enough for a hepatocyte to transform into a cholangiocytes, as then any hepatocyte in the tissue adjacent to a cholangiocytes would be able to transform at some point, causing a homogeneous spread of cholangiocytes, which is not observed experimentally. Hence, the second condition for the transformation is that the hepatocyte must delimit an existing (possibly small) lumen. Third, we also need to introduce a transition contact time T sig which specifies the minimum time the signalling needs in order for a hepatocyte to change to a cholangiocyte. A value T sig = 0 would mean here that the transformation is immediate, which is not observed as some hepatocytes may express weak SOX9+ signals but have a much larger size than a mature cholangiocyte. The signaling time is likely limited by the cell cycle time. When a hepatocyte gets the signal to transform, we assume that it will be fully differentiated upon the next cell division. We set here T sig = 2h as default value.

Initialization of model and boundary conditions of bile-duct system
The cells are initially positioned side by side on concentric circles with curvature determined by the radius of the portal vein. The most inner layer represents the portal vein endothelial for which we assume the positions are fixed during the simulations. The second layer are the mesenchyme, followed by the hepatoblasts, partially represented by the DCM and the outer layers fully represented by the CBM.
A background pressure P b of the CBM hepatoblasts needs to be exerted at the outer border of the segment as a result of radial tissue growth and resistance from nearby cells. This is achieved in the model by exerting an inward body force to the outermost border cells such that the inner hepatoblasts represented by the DCM have an average internal pressure of P b . For the DCM, the pressure is calculated by Eq.6. For the CBM the pressure p i on a cell i is derived from the virial stress and given by: being the stress tensor quantifying the stresses cell i experiences subject to contact forces F ij with other cells j [9]. Here, r ij is the vector pointing from the center of cell i to the cell j with || r ij || = d ij /2 and V i is the sampling volume which can be taken as the cell volume.
To ensure the cells remain in a planar configuration, we apply a small penalty force on the center-of-mass of all the cells each time the center of mass moves away from the XY plane. A simulation then starts in which only the cells positions are updated and a mechanical equilibrium is reached. This configuration is then used for all further simulations.

Calculation of the bile duct lumen area
We obtain the lumen area by considering the triangles of the cells that are marked by tracer particles. The procedure to compute the lumen volume is conceptually similar to how the volume of a cell is computed given its triangulated structure (see ref. [12]). First, the geometric center of all these triangles is computed, resulting in a vector pointing to the geometric center of the bile duct. From this point, we compute the signed volume of each tetrahedron with as base a marked triangle and as top the geometric center. The sign here is determined by the normal vector of the triangle. This is done for all marked triangles, and subsequently these volumes are then summed up. This gives a reasonable estimation of the volume confined by the marked cells.

Computational scheme executed during simulation
The simulations have been performed in TiSim, a computational platform for tissue simulations based on agent-based models. This software will be released in the near future.
In every timestep, the new positions of the nodes of the DCM and the centers of the CBM are sought, given the forces that act on them. To compute these forces, one must execute a contact detection algorithm over the arrays of interacting objects to know which objects (triangles, spheres) can interact. This is a computationally expensive task which requires an effective method. We have here used a bounding box algorithm associated which each triangle and sphere. To compute the cell positions, the system of equations (1) and (2) in the main text are solved simultaneously. This is a linear system which can be solved efficiently by a conjugate gradient (CG) method. As the CG iterations are taking a large part of the computational time, especially as the system becomes larger, we have parallelized some loops in the CG algorithm using OpenMP. The optimal number of threads is here typically 2 to 3. After the velocities are calculated, we use an Euler scheme to obtain the positions. We have used in all simulations a stable timestep of 1s. The full bile duct system of cells after 24h of simulation time contains about 20000 nodes from the DCM, 500 nodes from the CBM, and 7500 tracer particles. The time so solve this system takes about 48h on a modern machine (Intel(R) Xeon(R) Gold 6136).
Based on the cell positions and forces, actions are undertaken in the algorithm. In Fig A2 , a model flow chart is given with the principal algorithms that are executed every timestep. The algorithms that are cell -type dependent, are colored according to the cell type. Figure A2. Computational scheme executed during simulation.